
###########################################################################################
####### Main Analysis for "Legislative Effectiveness in the American States"
##########################################################################################
### Includes code to replicate all tables/figures in the paper, in order, EXCEPT:
### *** For Figure 1 and Table A4, see separate `SLES_analysis_nccppr.R` script 
### *** For all Appendix materials except Table A4, see separate `SLES_analysis_appendix.R` script 
###########################################################################################
### Originally run using R version 4.0.2 (2020-06-22)

#rm(list=ls())
options(stringsAsFactors = FALSE, scipen = 999)

### **********************************************************************
### Set Working Directory to Folder Containing Replication Files and Data
### **********************************************************************
working_directory <- ""

setwd(working_directory)

### Install/Load Packages
# install.packages(c("tidyverse", "estimatr", "ggridges", "glue", "stargazer"))
library(tidyverse) # version 1.3.0
library(estimatr)  # version 0.22.0
library(ggridges)  # version 0.5.2
library(glue)      # version 1.4.2
library(stargazer) # version 5.2.2

### Create Figures Directory in Working Directory if Not Present
if(!dir.exists("Figures")){
  dir.create("Figures")
}

### Create Tables Directory in Working Directory if Not Present
if(!dir.exists("Tables")){
  dir.create("Tables")
}

#######################################################
## Load Cleaned/Aggregated SLES Data and State Covariates
####################################################

#### Individual Level File - 1 row per legislator x legislative term
sles <- read_csv("Data/SLES_individual_level.csv", col_types = cols(district = 'c'))
glimpse(sles)

#### State Level File - 1 row per state x legislative term
sles_agg <- read_csv("Data/SLES_state_level.csv", col_types = cols())
glimpse(sles_agg)



#########################################################
### Figure 1: Criterion Validation in North Carolina
#########################################################


### See `SLES_analysis_nccppr.R` script



#########################################################
### Figure 2: Correlation between SLES vs Lagged SLES
#########################################################

## Pull Overall Correlations to Add to Plot: SLES T vs T-1 by Majority Status
sles %>% 
  group_by(sles_id, state, chamber) %>%
  arrange(term) %>%
  mutate(lagged_SLES = lag(SLES), 
         status = case_when(
           in_majority == 1 & lag(in_majority) == 0 ~ "Entered Majority",
           in_majority == 0 & lag(in_majority) == 1 ~ "Entered Minority",
           in_majority == lag(in_majority) ~ "No Change in Party Control",
           TRUE ~ NA_character_)) %>%
  filter(!is.na(status)) %>%
  group_by(status) %>%
  summarize(group_corr = cor(lagged_SLES, SLES),
            .groups = "drop")

### Compute Lagged Comparisons
fig_two_df <- sles %>% 
  group_by(sles_id, state, chamber) %>%
  arrange(term) %>%
  mutate(lagged_SLES = lag(SLES), 
         status = case_when(
           in_majority == 1 & lag(in_majority) == 0 ~ "Entered Majority",
           in_majority == 0 & lag(in_majority) == 1 ~ "Entered Minority",
           in_majority == lag(in_majority) ~ "No Change in Party Control",
           TRUE ~ NA_character_)
         ) %>%
  filter(!is.na(status))

#### Figure 2
figure_two <- fig_two_df %>%
  ggplot(aes(x = lagged_SLES, y = SLES, color = status, group = status)) + 
  geom_point(alpha = .15, size = .25) + 
  geom_smooth(method = "lm", formula = "y ~ x") + #  + I(x^2)
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") + 
  facet_wrap( ~ status) + 
  xlim(0, 10) + 
  ylim(0, 10) + 
  theme(legend.title = element_blank(), legend.position = "bottom") + 
  ylab("SLES (Current Term)") + 
  xlab("SLES (Prior Term)") + 
  ### Add Correlations as Text
  geom_text(aes(x, y, label=text), size = 4, vjust = 1, hjust = 0,
            show.legend = FALSE,
            data= data.frame(x=0.1, y=9.9,
                             text = c("\u03c1 = 0.38","\u03c1 = 0.40","\u03c1 = 0.66"), 
                             status = c("Entered Majority", "Entered Minority", "No Change in Party Control")))

figure_two

####### Save
ggsave("Figures/Figure2.jpeg", figure_two, height = 4.5, width = 9)

### To save as PDF with correct symbols:
# cairo_pdf("Figures/Figure2.pdf", height = 4.5, width = 9)
# print(figure_two)
# dev.off()


### Specifications to Replicate Results for Lines of Best Fit in Figure 2: 
summary(lm(SLES ~ lagged_SLES, data = fig_two_df %>% filter(status == "Entered Majority")))
summary(lm(SLES ~ lagged_SLES, data = fig_two_df %>% filter(status == "Entered Minority")))
summary(lm(SLES ~ lagged_SLES, data = fig_two_df %>% filter(status == "No Change in Party Control")))


rm(figure_two, fig_two_df)


#########################################################
### TABLE 1: Determinants of State Legislative Effectiveness Scores
#########################################################
# Follows Volden and Wiseman's 2014 LES Model + State/Term FEs
# ~ Term FEs done by Biennial Start Years. So any term that starts in 1995 or 1996 coded as 1995

### Models 1.1, 1.2, 1.3
base_model <- lm(SLES ~ seniority + comm_chair + in_majority + Leader_majleader + Leader_minleader + Leader_speakerpres + power_comm + ideo_med_distance + female + I(pred.race == "black") + I(pred.race == "hispanic") + vote_share + I(vote_share^2) + factor(state_chamber) + factor(biennial_grps), data = sles)
base_model_H <- lm(SLES ~ seniority + comm_chair + in_majority + Leader_majleader + Leader_minleader + Leader_speakerpres + power_comm + ideo_med_distance + female + I(pred.race == "black") + I(pred.race == "hispanic") + vote_share + I(vote_share^2) + factor(state) + factor(biennial_grps), data = sles, subset = sles$chamber == "lower")
base_model_S <- lm(SLES ~ seniority + comm_chair + in_majority + Leader_majleader + Leader_minleader + Leader_speakerpres + power_comm + ideo_med_distance + female + I(pred.race == "black") + I(pred.race == "hispanic") + vote_share + I(vote_share^2) + factor(state) + factor(biennial_grps), data = sles, subset = sles$chamber == "upper")

### Clustering SE's for use with Stargazer
base_SEs <- list(starprep(base_model, clusters = sles[-base_model$na.action,]$sles_id, se_type = "stata")[[1]], 
                 starprep(base_model_H, clusters = sles[sles$chamber == "lower",][-base_model_H$na.action,]$sles_id, se_type = "stata")[[1]],
                 starprep(base_model_S, clusters = sles[sles$chamber == "upper",][-base_model_S$na.action,]$sles_id, se_type = "stata")[[1]])

### Table 1 - Determinants of SLES
stargazer::stargazer(base_model, base_model_H, base_model_S, se = base_SEs,
                     type = "text", 
                     #type = "html", out = "Tables/Table1.doc",
                     dep.var.caption = "Dependent variable: SLES",
                     dep.var.labels.include = FALSE,
                     column.labels = c("Full Sample (1.1)",
                                       "Lower Chambers (1.2)",
                                       "Upper Chambers (1.3)"),
                     model.numbers = FALSE,
                     omit = c("state", "biennial"),
                     omit.stat = c("f", 'ser', 'adj.rsq'),
                     order = c("seniority", "chair", "in_majority", "majleader", "minleader", "speakerpres", "in_maj"),
                     covariate.labels = c('Seniority', "Committee Chair", "Majority Party", "Majority Leadership", "Minority Leadership", 
                                          "Speaker/President", "Power Committee", "Distance from Median", "Female",
                                          "African-American", "Hispanic", "Vote Share", "Vote Share\\^2", "Constant"),
                     add.lines = list(c("State-Chamber FE", "Yes", "No", "No"), 
                                      c("State FE", "No", "Yes", "Yes"),
                                      c("Term FE", rep("Yes", 3))),
                     star.char = c("+", "*", "**"),
                     notes = "+p<0.1; *p<0.05; **p<0.01", notes.append = FALSE)

rm(base_model, base_model_H, base_model_S, base_SEs)


#########################################################
#### Figure 3: Seniority by Majority/Minority/Chair
#########################################################

seniority_trends <- sles %>% 
  ### Filtering to Majority Party Non-Chairs and Minority Party
  filter(!(in_majority == 1 & comm_chair == 1)) %>%
  group_by(seniority, in_majority) %>% 
  summarize(median_SLES = median(SLES), .groups = "drop") %>%
  mutate(trend_group = case_when(in_majority == 1 ~ "Majority Party Non-Chair", in_majority == 0 ~ "Minority Party", TRUE ~ NA_character_)) %>%
  select(-in_majority) %>%
  bind_rows(
    ### Adding in Majority Party Chairs
    sles %>% 
      filter(comm_chair == 1 & in_majority == 1) %>% 
      group_by(seniority) %>% 
      summarize(median_SLES = mean(SLES), .groups = "drop") %>%
      mutate(trend_group = "Majority Party Chair")
  ) %>%
  filter(seniority <= 10)

figure_three <- 
  ggplot(seniority_trends, aes(x = seniority, y = median_SLES, color = trend_group, fill = trend_group)) +
  geom_point() + 
  geom_smooth(se = FALSE, method = "loess", formula = "y ~ x") + 
  ylim(0, 2) + 
  scale_x_continuous(breaks = seq(1, 11, 2)) + 
  xlab("Seniority (Terms in State Legislature)") + 
  ylab("Average SLES") +
  theme(legend.title = element_blank())

figure_three

####### Save
ggsave("Figures/Figure3.pdf", figure_three, height = 4, width = 8)

### Specifications to Replicate Loess Fit in Figure 3: 
predict(loess(median_SLES ~ seniority, data= seniority_trends %>% filter(trend_group == "Majority Party Chair")), newdata = 1:10)
predict(loess(median_SLES ~ seniority, data= seniority_trends %>% filter(trend_group == "Majority Party Non-Chair")), newdata = 1:10)
predict(loess(median_SLES ~ seniority, data= seniority_trends %>% filter(trend_group == "Minority Party")), newdata = 1:10)
# Below contains the exact x/y pairs used for plotting the line from this fit by ggplot
ggplot_build(figure_three)$data[[2]] 

rm(figure_three, seniority_trends)


#########################################################
#### Figure 4: Distance from the Median by Majority/Minority
#########################################################

figure_four <- sles %>%
  mutate(plot_group = case_when(in_majority == 1 ~ "Majority Party", in_majority == 0 ~ "Minority Party", TRUE ~ NA_character_)) %>%
  ggplot(aes(x = ideo_med_distance, y = SLES, color = plot_group, group = plot_group)) + 
  geom_point(alpha = .15, size = .1) +
  geom_smooth(method = "lm", formula = "y ~ x") + 
  ylim(0, 5) + 
  xlab("Distance from Chamber Median") + 
  ylab("Average SLES") +
  theme(legend.title = element_blank())

figure_four

#### Save!
ggsave("Figures/Figure4.pdf", figure_four, height = 4.5, width = 7)

### Specifications to Replicate Results for Lines of Best Fit in Figure 4: 
summary(lm(SLES ~ ideo_med_distance, data = sles %>% filter(in_majority == 1)))
summary(lm(SLES ~ ideo_med_distance, data = sles %>% filter(in_majority == 0)))


rm(figure_four)


##############################################################
#### Figure 5: Gender Differences in Effectiveness Scores Across States
################################################################


figure_five <- sles %>%
  mutate(st_factor = fct_rev(as.factor(state))) %>% 
  ### Order by Mean Cross-Chamber Effectiveness Difference (Male - Female)
  group_by(state, chamber) %>% 
  mutate(sort_var = mean(SLES[female == 0]) - mean(SLES[female == 1])) %>% 
  ungroup() %>% 
  mutate(st_factor = fct_rev(fct_reorder(as.factor(state), desc(sort_var), mean))) %>%
  mutate(gender = ifelse(female == 1, "Female", "Male"),
         chamber = ifelse(chamber == "lower", "Lower Chambers", "Upper Chambers")) %>%
  ggplot(aes(y = st_factor)) +
  geom_density_ridges(aes(x = SLES, fill = gender), alpha = .65, scale = 3, color = "white", from = 0, to = 4) +
  xlab("SLES") + 
  ylab("State") + 
  scale_fill_cyclical(
    breaks = c("Male", "Female"),
    values = c("yellow3", "dodgerblue"), # c("gray10", "gray80"),
    name = "", 
    guide = "legend") +
  theme_ridges(grid = FALSE) + 
  facet_wrap(~chamber) + 
  theme(axis.title.x = element_text(hjust = 0.5), 
        axis.title.y = element_text(hjust = 0.5), 
        legend.position = "right")


figure_five

#### Save!
ggsave("Figures/Figure5.pdf", figure_five, height = 9, width = 8)


### Exact densities used for plotting Figure 5 in ggplot
ggplot_build(figure_five)$data[[1]] 


rm(figure_five)


##############################################################
#### Figure 6: Comparing Baseline Relationships in Lower/Upper Chambers
################################################################

compare_dat <- sles %>% mutate(lower_chamber = ifelse(chamber == "lower", 1, 0))
compare_mod <- lm_robust(SLES ~ seniority*lower_chamber + comm_chair*lower_chamber + in_majority*lower_chamber + Leader_majleader*lower_chamber + Leader_minleader*lower_chamber + Leader_speakerpres*lower_chamber + power_comm*lower_chamber + ideo_med_distance*lower_chamber + female*lower_chamber + I(pred.race == "black")*lower_chamber + I(pred.race == "hispanic")*lower_chamber + vote_share*lower_chamber + I(vote_share^2)*lower_chamber + factor(state) + factor(biennial_grps),
                         data = compare_dat, 
                         clusters = sles_id, 
                         se_type = "stata")

chamber_interactions <- tidy(compare_mod) %>% 
  filter(grepl("lower_chamber:|:lower_chamber", term)) %>%
  mutate(term_fct = dplyr::recode(term,
                                  `seniority:lower_chamber` = "Seniority",
                                  `lower_chamber:comm_chair` = "Committee Chair",
                                  `lower_chamber:in_majority` = "Majority Party",
                                  `lower_chamber:Leader_majleader` = "Majority Leadership",
                                  `lower_chamber:Leader_minleader` = "Minority Leadership",
                                  `lower_chamber:Leader_speakerpres` = "Speaker/President",
                                  `lower_chamber:power_comm` = "Power Committee",
                                  `lower_chamber:ideo_med_distance` = "Distance from Median",
                                  `lower_chamber:female` = "Female",
                                  `lower_chamber:I(pred.race == "black")TRUE` = "African American",
                                  `lower_chamber:I(pred.race == "hispanic")TRUE` = "Hispanic",
                                  `lower_chamber:vote_share` = "Vote Share",
                                  `lower_chamber:I(vote_share^2)` = "Vote Share Squared")) %>%
  mutate(term_fct = factor(term_fct, levels = .$term_fct)) %>%
  mutate(term_fct = forcats::fct_rev(term_fct))

figure_six <- chamber_interactions %>%
  filter(!grepl("Vote", term_fct)) %>%
  ggplot(aes(x = estimate, y = term_fct)) +
  geom_point(size = .75) +
  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), colour="black", size = .35, height=.2) +
  geom_errorbarh(aes(xmin=estimate - qnorm(.95) * std.error, xmax = estimate + qnorm(.95) * std.error), colour="black", size = .65, height=0) +
  xlab("Difference in Point Estimates (Lower Chambers - Upper Chambers)") +
  ylab("") +
  geom_vline(xintercept = 0, linetype = "dashed")

figure_six

### Save
ggsave("Figures/Figure6.pdf", figure_six, height = 4, width = 7)

rm(compare_dat, compare_mod, chamber_interactions, figure_six)



##############################################################
#### Figure 7: Majority vs Minority Party Distributions by State
################################################################

figure_seven <- sles %>%
  mutate(st_factor = fct_rev(as.factor(state))) %>% 
  mutate(chamber = case_when(chamber == "upper" ~ "Upper Chambers",
                             chamber == "lower" ~ "Lower Chambers")) %>%
  ### Order by Mean Cross-Chamber Difference In Majority/Minority Party Effectiveness
  group_by(state, chamber) %>% 
  mutate(sort_var = mean(SLES[in_majority == 1]) - mean(SLES[in_majority == 0])) %>% ungroup() %>% 
  mutate(st_factor = fct_rev(fct_reorder(as.factor(state), desc(sort_var), mean))) %>%
  ### Drop NE -- Nonpartisan
  filter(state != "NE") %>%
  mutate(in_majority = ifelse(in_majority == 1, "Majority Party", "Minority Party"),
         chamber = tools::toTitleCase(chamber)) %>%
  ggplot(aes(y = st_factor)) +
  geom_density_ridges(aes(x = SLES, fill = in_majority), alpha = .65, scale = 3, color = "white", from = 0, to = 4) +
  xlab("SLES") + 
  ylab("State") + 
  scale_fill_cyclical(
    breaks = c("Majority Party", "Minority Party"),
    values = c("dodgerblue", "yellow3"), # c("gray10", "gray80"),
    name = "", 
    guide = "legend") +
  theme_ridges(grid = FALSE) + 
  facet_wrap(~chamber) + 
  theme(axis.title.x = element_text(hjust = 0.5), 
        axis.title.y = element_text(hjust = 0.5), 
        legend.position = "right")

figure_seven

## Save
ggsave('Figures/Figure7.pdf', figure_seven, height = 9, width = 8)

### Exact densities used for plotting Figure 7 in ggplot
ggplot_build(figure_seven)$data[[1]] 

rm(figure_seven)


#########################################################
### TABLE 2: Determinants of Majority-Party Advantage
#########################################################

### Models 2.1, 2.2
mm_diff_mod <- lm(maj_min_SLES_meddiff ~ SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + squire_index + aj_anymaj_setscal + aj_comm_gatekeep + unified_gov + log(chamber_size) + term_limits_enacted, data = sles_agg)
mm_overlap_mod <- lm(maj_above_min_SLES_median ~ SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + squire_index + aj_anymaj_setscal + aj_comm_gatekeep + unified_gov + log(chamber_size) + term_limits_enacted, data = sles_agg)

mm_SEs <- list(
  starprep(mm_diff_mod, clusters = sles_agg[-mm_diff_mod$na.action,]$state_chamber, se_type = "stata")[[1]],
  starprep(mm_overlap_mod, clusters = sles_agg[-mm_overlap_mod$na.action,]$state_chamber, se_type = "stata")[[1]]
)


### Table 2
stargazer(mm_diff_mod, mm_overlap_mod, 
          se = mm_SEs, 
          omit.stat = c("f", "ser", "adj.rsq"),
          type = "text",
          #type = "html", out = "Tables/Table2.doc",
          column.labels = c("SLES Partisan Difference (2.1)", 
                            "Share More Effective (2.2)"), 
          model.numbers = FALSE,
          dep.var.labels.include = FALSE,
          order = c("SM_", "partisan_seatshare", "squire", "^aj_", "^maj_"),
          covariate.labels = c("Polarization", "Majority Party Heterogeneity", "Minority Party Heterogeneity",
                               "Partisan Seat Share Imbalance", "Legislative Professionalism",
                               "Committee Gatekeeping Power", "Majority Party Controls Calendar", 
                               "Unified Government",
                               "Log Chamber Size", "Term Limits"),
          star.char = c("+", "*", "**"),
          notes = "+p<0.1; *p<0.05; **p<0.01", notes.append = FALSE)


rm(mm_diff_mod, mm_overlap_mod, mm_SEs)


#########################################################
### TABLE 3: Institutional Design and Patterns of State Legislative Effectiveness
#########################################################

### Models 3.1 - 3.4
mod_majmin <- lm(maj_min_SLES_meddiff ~ log(leg_salary_yr_avg + 1) + log(slength_avg_all) + ncsl_num_staff_per_leg + aj_anymaj_setscal + aj_comm_gatekeep + aj_chamb_vote_comm_appts + num_chamber_comms + log(chamber_size) + term_limits_enacted + SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + unified_gov, data = sles_agg) 
mod_chair <- lm(chair_SLES_meddiff ~ log(leg_salary_yr_avg + 1) + log(slength_avg_all) + ncsl_num_staff_per_leg + aj_anymaj_setscal + aj_comm_gatekeep + aj_chamb_vote_comm_appts + num_chamber_comms + log(chamber_size) + term_limits_enacted + SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + unified_gov, data = sles_agg) 
mod_majsenior <- lm(majfrosh_SLES_meddiff ~ log(leg_salary_yr_avg + 1) + log(slength_avg_all) + ncsl_num_staff_per_leg + aj_anymaj_setscal + aj_comm_gatekeep + aj_chamb_vote_comm_appts + num_chamber_comms + log(chamber_size) + term_limits_enacted + SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + unified_gov, data = sles_agg) 
mod_minsenior <- lm(minfrosh_SLES_meddiff ~ log(leg_salary_yr_avg + 1) + log(slength_avg_all) + ncsl_num_staff_per_leg + aj_anymaj_setscal + aj_comm_gatekeep + aj_chamb_vote_comm_appts + num_chamber_comms + log(chamber_size) + term_limits_enacted + SM_median_diff + SM_maj_sd + SM_min_sd + partisan_seatshare_diff + unified_gov, data = sles_agg) 

mm_SEs <- list(
  starprep(mod_majmin, clusters = sles_agg[-mod_majmin$na.action,]$state_chamber, se_type = "stata")[[1]],
  starprep(mod_chair, clusters = sles_agg[-mod_chair$na.action,]$state_chamber, se_type = "stata")[[1]],
  starprep(mod_majsenior, clusters = sles_agg[-mod_majsenior$na.action,]$state_chamber, se_type = "stata")[[1]],
  starprep(mod_minsenior, clusters = sles_agg[-mod_minsenior$na.action,]$state_chamber, se_type = "stata")[[1]]
)

### Table 3
stargazer(mod_majmin, mod_chair, mod_majsenior, mod_minsenior,
          se = mm_SEs, 
          omit.stat = c("f", "ser"),
          type = "text",
          #type = "html", out = "Tables/Table3.doc",
          column.labels = c("SLES Partisan Difference (3.1)", 
                            "SLES Chair Difference (3.2)", 
                            "Majority SLES Seniority Difference (3.3)", 
                            "Minority SLES Seniority Difference (3.4)"),
          model.numbers = FALSE,
          dep.var.caption = "Dependent variable: ",
          dep.var.labels.include = FALSE,
          covariate.labels = c("Log Annual Salary", "Log Session Length", "Staff per Legislator",
                               "Majority Party Controls Calendar", "Committee Gatekeeping Power", "Chamber Votes on Committee Appointments",
                               "Number of Committees", "Log Chamber Size", "Term Limits",
                               "Polarization", "Majority Party Heterogeneity", "Minority Party Heterogeneity",
                               "Partisan Seat Share Imbalance", "Unified Government", 
                               "Constant"),
          star.char = c("+", "*", "**"),
          notes = "+p<0.1; *p<0.05; **p<0.01", notes.append = FALSE)

rm(mod_majmin, mod_chair, mod_majsenior, mod_minsenior, mm_SEs)



